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Abstract. QCD at finite temperature and density is becoming increasingly im- 
portant for various experimental programmes, ranging from heavy ion physics to 
astro-particle physics. The non-perturbative nature of non-abelian quantum field 
theories at finite temperature leaves lattice QCD as the only tool by which we 
may hope to come to reliable predictions from first principles. This requires care- 
ful extrapolations to the thermodynamic, chiral and continuum limits in order 
to eliminate systematic effects introduced by the discretization procedure. After 
an introduction to lattice QCD at finite temperature and density, its possibili- 
ties and current systematic limitations, a review of present numerical results is 
given. In particular, plasma properties such as the equation of state, screening 
masses, static quark free energies and spectral functions are discussed, as well as 
the critical temperature and the QCD phase structure at zero and finite density. 

1 Introduction 

Quantumchromodynamics (QCD) is the theory of the strong interactions, describing nuclear 
matter in terms of its constituent quarks and gluons and their interactions. One of its key 
features is asymptotic freedom, i.e. the fact that the coupling strength is decreasing as a function 
of momentum transfer of an interaction. Thus, while the theory is amenable to perturbation 
theory at large momenta, it is non-perturbative for energy scales < 1 GeV and lattice QCD is 
the only known method for first principles calculations in this regime. 

At a critical temperature Tc w 200 McV, QCD predicts a transition between the famil- 
iar confined hadron physics and a deconfined phase of quark gluon plasma (QGP). At the 
same temperature, the weakly broken chiral symmetry, responsible for the three light pions, 
gets restored. QCD at high temperatures is of outstanding importance in today's theoretical 
and experimental nuclear and particle physics programmes. A thermal environment with suf- 
ficiently high temperatures for a QCD plasma has certainly existed during the early stages of 
the universe, which passed through the quark hadron transition on its way to its present state. 
Current and future heavy ion collision experiments are recreating this primordial plasma at 
RHIC (BNL), LHC (CERN) and FAIR (GSI). These studies will have a bearing far beyond 
QCD in the context of early universe and astro-particle physics. Many prominent features of 
the observable universe, such as the baryon asymmetry or the seeding for structure formation, 
have been determined primordially in hot plasmas described by non-abelian gauge theories. The 
QCD plasma serves as a prototype for those, since it is the only one we can hope to produce 
in laboratory experiments. Moreover, certain QCD plasma properties entering calculations of 
dark matter abundances need to be known at percent level accuracy in order to match the ever 
more precise astro-physical data expected in the near future [1]. On the other hand, for high 
densities and low temperatures, exotic non-hadronic phases such as a color superconductor have 
been predicted [2], and such physics might be realized in the cores of compact stars. 
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For these applications we need to know how the properties of QCD change under extreme 
conditions of temperature and density. This entails a determination of the QCD phase diagram 
and the associated critical properties, a quantitative understanding of the equation of state, the 
way in which hadron properties get modified as well as the nature of the lightest excitations in 
non-hadronic phases. 

Early hopes that the high temperature plasma phase might be accessible to perturbation 
theory have proved wrong, as we shall see. Indeed, today one rather speaks of the strongly cou- 
pled quark gluon plasma (sQGP). This leaves lattice QCD as the main computational tool even 
in that regime. However, lattice simulations are still struggling with limitations and systematic 
errors. For example, only in the last few years has it become possible to perform simulations for 
QCD at small baryon density, while the case of cold and high density matter remains unacces- 
sible to date. These lectures are intended to provide an introduction to the problems tackled by 
lattice simulations, their potential and limitations, as well as an overview over current results. 



2 Finite temperature QCD in the continuum and on the lattice 

2.1 QCD at finite temperature and density 

The thermodynamics of QCD is most conveniently described by the grand canonical partition 
function [3] 

Z{V,fif,T;g,mf,) = Tr(c-'^"-''f'^f^^^^ = J ZJADV^D?/. e"^"!^^! e"^^!'^''^'^^!, (1) 
with the euclidean gauge and fermion actions 

l/T 

Sg[Af^] = J dxo j '^^^ ^Tr Ff^^Ff_,„, 

V 

5'/[V', -0,^^] / dxQ I Sx^i}f{-i^,D^ + mf - ^if^Q)^f. (2) 
V 

The partition function depends on the external macroscopic parameters T,V^fif, as well as on 
the microscopic parameters like masses and the coupling constant. The index / labels the differ- 
ent quark flavours, and the conserved quark numbers corresponding to the chemical potentials 
jif are Qf = ipjoi^- In the following we will consider mostly two and three flavours of quarks, 
and always take m„ = nid- The case — rriu^d is then denoted by Nf — 3, while Nf — 2 + 1 
implies ^ ruu^d- For our purposes we will couple all flavours to the same chemical potential 
/i unless otherwise stated. The chemical potential for baryon number is then /is = 3/i. Once the 
partition function is known, thermodynamic properties such as free energy, pressure, average 
particle numbers or the thermal expectation value of an operator O readily follow, 

F = -TlnZ, P=^i^-^, (Q^)^^iZ^, (0)=Z-iTr(Oe-(^-««/)/^). (3) 
uV ^M/ 

In the thermodynamic limit V — > cx), we are interested in the corresponding densities, i.e. f = 
F/V,p^P/V^~f,.... 

While the formulation of QCD thermodynamics is straightforward, controlled evaluations 
are prohibitively difficult. Perturbation theory, the tool so successful for electroweak physics 
at zero temperature, fails us completely in this context. The reasons are two-fold. Firstly, the 
running coupling is still not weak enough at temperature scales of interest, (/(lOO— 300MeV) ^ 1. 
Worse still is the fact that finite T perturbation theory for non-abelian gauge theories displays 
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an insurmountable infrared problem 4J. At finite T, a gauge boson self coupling gets dressed 
with a thermal distribution function which for small momenta behaves as 

-wf-, ' — • (4) 

Thus, the effective expansion parameter diverges for perturbatively zero mass particles such 
as gluons, no matter how weak the coupling, i.e. this problem also exists in the symmetric 
electroweak phase. Higher orders do generate an effective gluonic mass scale m g^T to cure 
this divergence, and resummation methods have been devised to calculate it ^5 . But now the 
coupling cancels out of the expansion parameter and all loop orders contribute in the same 
way, leaving convergence of such series unclear. Hence, a fully non-perturbative treatment of 
the problem is warranted. 



2.2 The Lattice Formulation 

The thermal quantum field theory is discretized by introducing a euclidean space-time lattice 
X Nt with lattice spacing a, such that volume and temperature are 

V^iaLf, T.-i-. (5) 

The fermion fields live on the lattice sites x, whereas the gauge fields are represented by link 
variables Ufj,{x) e S'[/(3) connecting the sites. The simplest actions used in many thermody- 
namics calculations are the Wilson gauge action and the staggered, or Kogut-Susskind, fermion 
action 

Sa[U]^Y. E /^fl-^ReTrf/^y Sf'^Y.^{x)A'Lyx{y). (6) 

respectively. Here Up = U ^{x)U,y{x + ajl)Uj^{x + av)Ul{x) is the elementary plaquette and the 
bare lattice and contimmm gauge couplings are related by (3 = &/g^. The staggered fermion 
matrix is given by 

1 ^ 

+ \{5.+o,yU,{x)e^^^^-5.^^^^^U,{y)^c-^n + S^yc^rnj . (7) 

Note that the chemical potential a/i enters through the temporal links only [6 . Now the Gaus- 
sian integral over the quark fields can be performed, leading to the partition function 

Z{L,afi,Nt;P,Nf,amf) ^ f DU YlidetMimj, e'^^^^l (8) 

The lattice action is not unique, any action reproducing the continuum action in the limit 
a ^ is in principle admissible. We refer to textbooks [7] for discussions of various actions and 
their advantages/disadvantages, as well as the daunting task of a chiral fermion formulation. 
While the above staggered fermions have a remnant U(l) chiral symmetry, they also exhibit 
spurious flavours (called tastes) which then have to be removed by taking the fourth root 
of the fermion determinant in Eq. Q. It is not yet settled whether this produces potentially 
hazardous non- localities or singularities, a subject of much current debate . Things empirically 
seem to work well as long as continuum limits are taken before chiral limits. On the other 
hand, Wilson fermions do not require roots of determinants, but they break chiral symmetry 
completely at finite lattice spacing and require additive quark mass renormalization. They 
also feature so-called exceptional configurations with negative eigenvalues of the determinant, 
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making simulations with light quarks very difficult. There are also fermion actions avoiding 
those issues, such as domain wall fermions or overlap fermions, for reviews see [9110) . However, 
at the present stage these actions are an order of magnitude more costly in simulation time and 
hence not yet widely used for dynamical thermodynamics simulations. It is for this reason that 
most results quoted in this review have been produced with staggered or Wilson fermions and 
their improved variants. The non-uniqueness of lattice actions can be exploited to construct 
actions which remove the 0(a), 0{a?) . . . corrections to continuum results, thus improving the 
approach to the continuum limit. Introductions to improvement can be found in [llj . 



2.3 Constraints on lattice simulations and systematic errors 

It is important to realize from the outset that current lattice simulations, and especially those 
at finite temperature and density, are still hampered by systematic errors and uncertainties. 
Let us discuss the origins of those. The Compton wave length of a hadron is proportional to its 
inverse mass, ~ J^i^^, and the largest of those constitutes the correlation length of the statistical 
system. To keep finite size as well as discretization errors under control, we need to require 

a < m^^ < aL. (9) 

Typical lattice sizes today are (12 — 32)'^ x 4, x 8 etc., depending on computing budgets and 
machines. For temperatures around T ~ 200 MeV, Eq. ([5]) with Nt — A — Q implies a ^ 0.1 — 0.3 
fm and aL ~ 1.5 — 3 fm. For the low T phase, Eq. ^ then tells us that the lightest hadron 
needs to be m^r > 250 MeV. The push to do physical quark masses is only just beginning to be 
a possibility on the most powerful machines and with the cheapest actions [i.e. staggered, with 
Wilson rapidly catching up). On the other hand, at high T screening masses scale as ran ~ T 
requiring N^"^ ^ 1 <C LN^^. Hence, the spatial lattice size should be significantly larger than 
the temporal one. This limits the feasible temperatures to several Tc- 

Once the simulations have been carried out, all quantities appear in lattice units, i.e. as 
dimensionless numbers, and need to be translated to physical units. This procedure of "setting 
the scale" introduces additional systematic errors, which are often much larger than the statis- 
tical errors of the simulations. Suppose we want to convert some measured critical temperature 
to physical units. It thus needs to be related to a quantity whose value in nature we know 
and which can also be obtained from a lattice simulation, for example a hadron mass at zero 
temperature, 

^ - ^ (10) 
TTiH NtianiH)' 

In practice this is difficult, because QCD with physical parameters is not quite doable numeri- 
cally yet. Moreover, one is often interested in the limit of heavy or zero quark masses, or with 
two or three degenerate flavours. In such cases one uses quantities related to the static potential, 
such as the string tension a or the Sommer scale rg |12j . which strictly speaking only exist in 
the pure gauge or heavy quark effective theories, 



T 1 ^ ro odV(r) 

425MeV; Tra — ^ with ^ ' 



[ay/^Nt' ' aNt dr 



1.65. (11) 



r=ro 



The values are provided through analyses of heavy quark effective theories applied to 6-quark 
systems, e.g. [13] . Fortunately the string tension turns out to be rather insensitive to the 
presence of light quarks, i.e. it is approximately the same independent of the quark content of 
QCD. 

Another difficulty arises for extrapolations to the continuum limit. Given some temperature 
T, the lattice spacing is set via T — l/{aNt), which then controls the gauge coupling (3 = 
6/g^{a). Thus, changing p changes temperature for a lattice with fixed A't. As a consequence, 
the cut-off in a simulation at fixed Nt varies as a function of temperature. E.g., if a simulation 
is performed for a set of quark masses fixed in lattice units, am/ = const, then changing 
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a(T) implies that the quark mass in physical units changes as well! In principle this can be 
avoided by simulating along the "lines of constant physics" , i. e. tuning bare masses together 
with the lattice spacing such that the masses in physical units stay constant. In practice, this 
is a formidable task since the lines of constant physics are, of course, not known and have to 
be mapped out non-perturbatively first. Again, simulations along lines of constant physics are 
only just being started. 



2.4 Finite baryon density and the sign problem 

As soon as a chemical potential for quark number is switched on, /i 7^ 0, a Monte Carlo 
evaluation of the partition function eq. ([8]) is essentially impossible due to the so-called sign 
problem of QCD. The problem is encapsulated in the 75-hermiticity of the Dirac operator, 

- 75^(-A^*)75, (12) 

which implies that detM is complex for the gauge group SU{3) and /i ^ 0. Note that the 
partition function and all physical observables remain real, i.e. the imaginary parts of the de- 
terminant cancel out once they are averaged over all gauge configurations. However, in a Monte 
Carlo evaluation using importance sampling, the determinant is evaluated in the background 
of single gauge configurations and interpreted as part of the probability weight for that config- 
uration. This is impossible if the determinant is complex. Methods to circumvent this problem 
will be discussed in sec. [5] 



2.5 The quenched limit, the chiral limit and physical QCD 

Much of our intuition regarding the quark hadron transition is based on the limiting cases of 
infinite or zero quark masses, i.e. the quenched (or static quark) and chiral limits, respectively. 
In these cases there exist true order parameters and the phase structure and transition can be 
discussed qualitatively based on symmetry breaking and universality arguments. 

For infinitely heavy quarks QCD reduces to Yang-Mills theory in the presence of static 
sources. The static quark propagator is given by a Wilson line closing through the time bound- 
ary, i.e. the Polyakov loop, 

Nt 

Li^)^Y[Uoix). (13) 

Xo 

The action displays a global symmetry under center transformations Uo{x) — > ZnUo{x), with 
Zn = expj27rn/3 S ^(3). On the other hand, static sources transform non-trivially, L(x) — > 
ZnL{yL). The free energy of a static quark is given by [2] 

(TrL(x)) - e~T^, (14) 

and constitutes an order parameter for confinement, i.e. (L) — for T < Tc, and (L) ^ for 
T > Tc- Hence there must be a true non-analytic phase transition corresponding to the breaking 
of the Z{?)) symmetry. One can construct an effective theory for the order parameter reflecting 
its Z{?)) symmetry |15l . and universality arguments suggest a first order phase transition. This 
is indeed borne out by simulations. 

In the limit of zero quark masses the classical QCD Lagrangian is invariant under global 
chiral symmetry transformations, the total symmetry being C/^(l) x SULiNf) x SUii{Nf). The 
axial J7yi(l) is anomalous, quantum corrections break it down to Z{Nf). The remainder gets 
spontaneously broken to the diagonal subgroup, SUL{Nf) x SUR{Nf) SUv{Nf), giving rise 
to iV| — 1 massless Goldstone bosons, the pions. The order parameter signalling chiral symmetry 
is the chiral condensate. 
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It is nonzero for T < Tc, when chiral symmetry is spontaneously broken, and zero for T > Tc- 
Effective models for the order parameter in this case are thus 0{2Nf) linear sigma models ;16j. 
The finite temperature phase transition then corresponds to chiral symmetry restoration. 

QCD with physical quark masses obviously does not correspond to either limit. The Z{3) 
symmetry is explicitly broken once there are dynamical fermions, and the chiral symmetry is 
broken by non-zero quark mass terms. Nevertheless, physical QCD displays confinement as well 
as three very light pions as "remnants" of those symmetries. In the presence of mass terms 
there is no true order parameter, i.e. the expectation values of the Polyakov loop as well as the 
chiral condensate are non-zero everywhere. Hence the deconfined or chirally symmetric phase 
is analytically connected with the confined or chirally broken phase, and there is no need for 
a non-analytic phase transition. The following questions then arise, which should be answered 
by numerical simulations: for which parameter values of QCD is there a true phase transition, 
and what is its order? Are confinement and chirality changing across the same single transition 
or are there different ones? If there is just one transition, which is the driving mechanism? If 
there is only a smooth crossover, how do the properties of matter change? 



2.6 Effective high T theory: dimensional reduction 

As we have seen, systematics presently constrain simulations of lattice QCD to temperatures 
below a few times T^. However, for applications to early universe physics, as well as our own 
understanding and comparisons with perturbation theory, we would also like to have non- 
perturbative information at much higher temperatures. This can be achieved by means of 
effective field theories. At large T, when the gauge coupling g{T) is sufficiently small, a hierarchy 
between different relevant scales of thermal QCD develops, 

27rT > > g^T. (16) 

The lowest non-vanishing bosonic Matsubara mode ~ 27tT is characteristic for non-interacting 
particles. The dynamics generates the Debye scale itie ^ gT, over which color electric fields 
are screened, and its non-perturbative analogue tum ^ g'^T for color magnetic fields Ai [4]. 

For physics on scales larger than the inverse temperature, |x| ~ ^/gT ^ this al- 

lows an effective theory approach in which the calculations are factorized: integration over the 
hard modes may be performed perturbatively by expanding in powers of the ratio of scales 
gT / (2'itT) ~ g/{2'K). This includes all non-zero Matsubara modes, in particular the fermions. 
It results in a 3d effective theory for modes ~ gT and softer, 

Seff= I dPx I i Tt{F,jF,j) + Tr(A^o)' + rn| Tr{Al) + \:,{TT{Alf^ , (17) 

and is known as dimensional reduction |17| . Since 4d euclidean time has been integrated over, 
Aq now appears as an adjoint scalar, and the effective parameters arc functions of the original 
ones, gl = g-'T, mE{N, N f, g, m{) ^ gT,\:i{N,Nf,g,m{)^ g^T. 

Discretization and simulation of this reduced problem is easy. Without fermions and in 
one dimension less, much larger volumes and finer lattices can be employed. Moreover, 3d 
gauge theories are superrenormalizable and the coupling scales linearly with lattice spacing. 
Hence, very accurate continuum limits can be obtained and systematic errors from simulations 
can be eliminated. However, the reduction step entails two approximations: the perturbative 
computation is limited to a finite order in g and neglects higher-dimensional operators, which 
are suppressed by powers of the scale ratio. The reduction step has been performed up to two- 
loop order [18j at which parameters have relative accuracy ©(g'*), while for correlation functions 
the error is [19] 5C/C ^ 0{g^). In the treatment of the electroweak phase transition, this is 
less than 5%[30], for QCD it depends on the size of the coupling g{T) and thus on T. It has 
been non-perturbatively checked by comparing simulations in 4d with those of the 3d effective 
theory, that the latter accurately reproduces static correlation lengths at temperatures as low 
as T > 2Tc |21I22| , thus allowing for a straightforward treatment of very large temperatures as 
well as detailed dynamical investigations in the plasma phase. 
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Fig. 1. The order parameters for deconfinement (L) (left) and chiral symmetry restoration (f/'^) 
(middle), together with their susceptibilities as a function of the gauge coupling for two-flavour QCD. 
From 23 . Right: Quark mass dependence of Tc for Nf — 2,3 with improved (p4) and standard 
staggered quarks. From [28] . 



3 Numerical results at finite T and /i = 

3.1 The (pseudo-) critical temperature 

The first task when investigating finite temperature QCD is to find the phase boundary. Thus, 
for QCD with a fixed quark content, we are interested in the critical (or pseudo-critical) tem- 
perature where the transition from the confining regime to the plasma regime takes place. 
The method to locate a transition usually is to look for rapid changes of observables like the 
Polyakov loop L, the chiral condensate or the plaquette, and for peaks of their fiuctuations, 
i.e. the generalized susceptibilities 

Xo = VNt{{0-{0})^). (18) 

The locations of these peaks define (pseudo-)critical couplings Pc which can be turned into tem- 
peratures through the knowledge of a zero temperature quantity like a hadron mass amniPc) 
at the same coupling, Tc/niH — {NtamHiPc))- In practice, often the two-loop beta function 
is used as a short cut, although this becomes valid only when the lattice spacing is fine enough 
to be in the perturbative regime. 

On coarse lattices Nt = 4, 6, all numerical evidence is consistent with confining and chiral 
properties changing at approximately the same (3c. An example is shown in fig. [T] (left and 
middle) for two-flavour QCD. On finite volumes, true non- analytic phase transitions do not 
exist and the identified phase boundary represents only a crossover. The nature and critical 
properties of the transition can be obtained from scaling behaviour of various quantities with 
volume and quark mass, as will be discussed in sec.[4l 

Tc for the pure SU{?>) gauge theory has been known rather precisely for quite a while [24l25j . 
It has since then been confirmed by studies using various improved actions, thus reducing finite 
lattice spacing effects. The number is most readily given in terms of the string tension, 

Tc/Vo^ = 0.632 ± 0.002, (19) 

which was obtained as a weighted average over all available lattice data [26]. For QCD with 
dynamical quarks, one finds Tc to be increasing with quark mass, as shown in fig. [T] (right). The 
data correspond to one lattice spacing only, i.e. are not yet contimmm results. Extrapolating 
the results to the chiral limit, one obtains 

9-fl CiCVt T-/(^'^^=^'^) MeV, clover-improved Wilson [57] 
- — navoi LjUU . - I (^73 _^ jyjgy^ improved staggered [21] 

3 - flavor QCD : = (154 ± 8) MeV, improved staggered [55] 
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Fig. 2. Left; Tc in units of ro as a function of pseudo-scalar mass for lattices with Nt — 4,6. The 
vertical line shows the physical value mpsro — 0.321(5). From [29]. Right: Continuum extrapolations 
of Tc, determined from susceptibilities of the chiral condensate, the strange quark susceptibility and 
the Polyakov loop for physical quark masses. From [3U] . 
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where mp has been used to set the scale. In view of the systematic errors introduced by the 
discretization of the fermion sector, it is reassuring that consistent results are obtained from 
staggered and Wilson fermions. Qualitatively we see that adding light degrees of freedom reduces 
the critical temperature for the transition. 

Having studied the flavour and quark mass dependence of the phenomenologically im- 
portant quantity is of course Tc for physical quark masses in the continuum limit. Two recent 
works have pushed in this direction |29|30j . and the comparison of the two is an interesting 
illustration of systematic effects in lattice calculations. The authors of 29] have worked with 
two lattice spacings Nt = 4,6, with the strange quark mass tuned to its physical value and a 
range of light quark masses, arriu^d = 0.1 — 0.4 TOs. Critical couplings are determined from peaks 
of chiral and Polyakov loop susceptibilities and found to be consistent with each other. They are 
then converted to temperature via a renormalization group improved two-loop beta function 
and the scale is set by tq (cf. sec. 12. 3|) . Finally, a combined extrapolation to the continuum and 
chiral limits is performed, as shown in fig. [2] (left), leading to the prediction Tc = 192(7)(4) 
MeV. 

A different number is quoted in |30j. In this work four lattice spacings are considered, 
Nt — 4, 6, 8, 10, for which the quark masses are tuned along the lines of constant physics, 
such that the meson masses correspond to nearly their physical values on all lattices. The 
critical couplings are determined through the chiral and strange quark susceptibility, as well 
as from the steepest change of the Polyakov loop. This is done for all lattice spacings, followed 
by a continuum extrapolation, as shown in fig. [2] (right). The scale in this work is set by 
the kaon decay constant /k, and the final value is based on the extrapolation of the chiral 
condensate, Tc — 151(3)(3). There are two striking observations to be made. Firstly, the Nt — 4 
data do not fall onto the straight a^-extrapolation line to the continuum limit, i.e. are not in 
the scaling region yet. Second, when finer lattices are considered, different observables lead 
to different critical temperatures. This is to be expected for a crossover, where Tc is only 
pseudo-critical and not uniquely defined. As we shall see in sec. lU this is the situation for 
physical QCD. Nevertheless, the gap is larger than expected, and also the ordering is puzzling: 
it suggests that there is a temperature range still displaying confinement while chiral symmetry 
is already restored, which goes against conventional wisdom. As a potential source for the 
discrepancy with [29], the authors of [30] identify cut-off effects, as is illustrated in fig. [3l If 
the continuum extrapolation is performed with two different methods to set the scale, then the 
data of [29^ appear to lead to inconsistent results. Likewise, a continuum extrapolation using 
only iV( = 4, 6 from the data in [30 would give a larger Tc than using the full data set. On the 
other hand, comparing the results obtained from the analysis of Polyakov loops between the two 
works gives only a small discrepancy. Clearly, cross checks with higher precision, more lattice 
spacings and other fermion discretizations are needed in order to come to fully conclusive results 
for the continuum. This discussion elucidates the difficulty of extracting accurate continuum 
information from lattice data, when statistical errors only appear to be small and under control. 
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Fig. 3. Left: Reanalysis of the data from [29] with different ways of setting the scale. Right: Same for 
the data from [30]. Both plots from [30j . 



3.2 The equation of state 

Energy density e{T) and pressure p{T) as a function of temperature are certainly among the 
most fundamental thermodynamic quantities of QCD governing, e.g.^ the expansion of the 
plasma in the early universe as well as in heavy ion collisions. The partition function for gases 
of free fermions and bosons are known exactly, and we recall the results in the relativistic and 
non-relativisic limits, the former corresponding to the famous Stefan-Boltzmann result: 

Relativistic Boson, m x (Fermion ) Non-relativistic, m ^ T 

Pr^9^T^ (i) p„. = ffr(^)^exp(-m/r) (20) 

Pr 9"^^ (s) Pjir "q^Prir Pnr 

For the fully interacting case, we need to compute the free energy density / = —T/VlnZ{V,T) 
of the QCD partition function, and the pressure follows as p = — /. A technical obstacle here 
is that, in a Monte Carlo simulation, one cannot compute the partition function directly, since 
that provides the probability weights and is normalized to one. What one can compute are 
expectation values (O). The most frequently used detour is called the integral method f3Tj, 
which computes a suitable derivative of the free energy resulting in a non-trivial expectation 
value, and then integrates the result. 

Note that this introduces a lower integration constant, which needs to be fixed for the result to 
be meaningful. While we do not know /(/3o) from first principles, we can choose (3q corresponding 
to a temperature below the phase transition, where the free energy should be well modelled by 
a weakly interacting glueball gas. Glueballs are heavy (> 1.6 GeV), and according to eq. ([20]) 
they produce exponentially small pressure which can be approximated by zero. 

Another difficulty is that strong discretization effects are to be expected. At high temper- 
ature the relevant partonic degrees of freedom have momenta of order ttT ~ 7r/(aiVf), i.e. on 
the scale of the lattice spacing, which thus strongly affects them. For the equation of state it is 
therefore particularly important to gain control over these effects and carry out the continuum 
limit a — > 0. This motivates the use of improved actions, designed to minimize cut-off effects in 
the approach to the continuum. 

The results of a computation of the pressure with an improved action [32j are shown in 
fig. lU The data have been obtained for Nf = 2,3 with (bare) mass nig/T = 0.4 as well as for 
Nf = 2-1-1 with a heavier mass rri'^/T = 1. For comparison the figure includes the continuum 
extrapolated pure gauge result. The figure shows a rapid rise of the pressure in the transition 
region. The critical temperature as well as the magnitude oip/T'^ refiect the number of degrees 
of freedom liberated at the transition. This last conclusion is firm, since the pressure also rises 
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Fig. 4. Left: Flavour dependence of the pressure for A^'t = 4 lattices compared to a continuum ex- 
trapolated pure gauge result. From [3^. Middle: The pressure for N j =2 + 1 with physical quark 
masses on A^t = 4,6. From jiSj. Right: Various orders of perturbation theory and computation within 
the 3d effective field theory to order ln(l/(7). A missing contribution to Oicf") has been adjusted by 
matching to the 4d lattice results, represented by the diamonds. From [35j . 



for fixed temperature when light quarks are added to the theory, consistent with the behaviour 
in the Stefan-Boltzmann limit. 

An important question then is whether these features survive in the continuum limit. Again, 
in pure gauge theory this can be firmly established by numerical extrapolation. First steps 
towards a continuum extrapolation of dynamical simulations have been taken in |33|34j . with 
consistent results. An example is shown in fig. [3] (middle), which was computed with bare quark 
masses tuned to their physical values. The trend to fall short of the Stefan-Boltzmann limit 
turns out to be confirmed on a finer lattice, and hence appears to be a robust result about the 
quark gluon plasma at moderate temperatures. 

Investigations of the equation of state thus provide us with an intriguing picture of the com- 
plexity of the quark gluon plasma. The sudden rise in pressure across the critical temperature 
has to be interpreted as a release of degrees of freedom, and in this sense it is justified to speak 
of a deconfining transition. On the other hand, the deconfined phase clearly shows remnant 
interactions and the released degrees of freedom are not appropriately characterized as a nearly 
free parton gas. 

At asymptotically high temperatures, the pressure eventually must meet the free gas limit 
because of asymptotic freedom. This can be studied in the framework of the dimensionally 
reduced effective high T theory. The pressure can be computed perturbatively to the order 

before the prohibitive Linde problem sets in at 0(5^) Fig. [4] (right) shows the poor 
convergence behaviour of the series up to that order. The non-perturbative contribution 
is dominated by the magnetic modes ~ g^T ■ However, it can be computed on the lattice 
by simulations of the 3d effective high temperature theory |36j . This requires matching of 
the coefficients between the effective and the full theory and converting from lattice to MS 
regularization schemes at four loop level, a task that has been recently completed However, 
there is one last missing contribution to 0((7^) coming from the perturbative scale ~ T. In the 
figure this coefficient has been chosen such that the results for the pressure computed in this 
framework match on to the lattice results at lower temperatures. 



3.3 Screening masses 

Essentially all static equilibrium properties of a thermal quantum field theory are encoded in 
its equal time correlation functions. These are quantities that are well defined and calculable 
to good precision by lattice methods. Unfortunately, these quantities are not directly accessible 
in heavy ion collision experiments. Nevertheless, their theoretical knowledge provides us with 
the relevant length scales in the plasma, from which conclusions about the active degrees of 
freedom and their dynamics may be drawn. The connected spatial correlation functions of 
gauge-invariant, local operators A{x)^ 

C(|x|) = (A(x)A(0)),^e-*^W, (22) 
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Fig. 5. Left: Screening masses for the pure gauge theory from A'^t = 4 lattices, corresponding to the 
continuum Oljl^ (circles) and Ol~ (squares) channels. From [41]. Right: Mesonic screening masses in the 
quenched chiral limit. Below Tc, Mp/Tc [42] is plotted, the line denotes the T — value. Open Wilson 
symbols denote anisotropic lattices [43]. Staggered pion data are extrapolated to a = [44], Wilson 
and staggered rho are from Nt = 8 |42) and 16 lattices. The free quark limit has not been corrected for 
finite volume effects. From 1261. 



fall off exponentially with distance. The "screening masses" M are the eigenvalues of the space- 
wise Hamiltonian of the corresponding lattice field theory, and classified by its symmetries. 
Because of the shortening of the euclidean time direction at T > 0, the continuum rotation 
symmetry of the hypertorus orthogonal to the correlation direction is broken down from 0(3) 
to 0(2) X Z{2), and its appropriate subgroup for the lattice theory is Df^. The irreducible rep- 
resentations and the classification of operators have been worked out for pure gauge theory 
[38 39J as well as for staggered quarks [40]. Physically, the screening masses correspond to the 
inverse length scale over which the equilibrated medium is sensitive to the insertion of a static 
source carrying the quantum numbers of A. Beyond 1 /M, the source is screened and the plasma 
appears undisturbed. Technically, the computation of these quantities is equivalent to spectrum 
calculations at zero temperature. 

Fig. [5] (left) shows results for the lowest lying screening masses corresponding to glue- 
ball operators around Tc- In the range 0.8 Tc < T < Tc, the lowest scalar screening mass 
is observed to be roughly 20% lower than the lightest scalar glueball at zero temperature, 
M{T)/Mc{T = 0) ~ 0.8. At Tc a sharp dip is observed, after which the screening masses 
appear to be proportional to T. Screening states with larger masses show the same qualitative 
behavior above Tc. 

Apart from pure glue operators A in eq. (|22p . also correlations of meson operators have 
been investigated, both in the quenched approximation as with various numbers of dynamical 
fcrmions. The picture which has emerged so far is illustrated by some of the available data shown 
in fig. [5] (right) . Below Tc, the screening masses do not show a marked temperature dependence 
and roughly agree with the corresponding zero temperature masses. At temperatures (slightly) 
above Tc spatial (as well as temporal) correlation functions reflect the restoration of the chiral 
SUL{Nf) X SUii{Nf) symmetry. In particular, the vector and axial vector channel become 
degenerate independent of the discretization and of the number of dynamical flavors being 
simulated. Moreover, the pion ceases to be a Goldstone boson and acquires a screening mass. 

At sufficiently high temperature, the screening masses are expected to approach the limit 
of freely propagating quarks, M 2ttT . In fact, already at temperatures of about 1.5 Tc, the 
results for the vector channel are not far from this value. If the finite volume effects for free 
quark propagation are taken into account, the deviations amount to about 15 % and decrease 
slowly with temperature. Thus, the hard fermionic modes ~ ttT behave quasi-perturbatively. 
Quenching effects are found to be below 5% for T > Tc [45], so the computer time saved on 
dynamical fermions can be invested in finer lattices and results are close to continuum physics. 
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Fig. 6. Pure gauge screening masses from simulations of the effective 3d theory in various quantum 
number channels for T — 2Tc (left) and T = lO^'^Tc (middle). Filled symbols refer to operators made 
of only magnetic gauge fields Ai ~ g^T, while open symbols also contain electric gauge fields Ao ~ gT. 
From 21 . Right: Interpolated temperature dependence of the two lowest Q'Y^ and the lowest OZ"*" 
states [46]. 



In the plasma phase, screening masses can also be investigated within the dimensionally 
reduced theory, cf. section 12.61 In this framework they correspond to the spectrum of the 
transfer matrix for the 3d theory. The associated Hamiltonian respects S0(2) planar rotations, 
two-dimensional parity P, -reflections R, and the symmetry is again 5*0(2) x Z{2) x Z{2) = 
0(2) X Z{2). Remember, however, that in this setup one is interested in soft modes, while the 
Matsubara frequency ^ 2'kT represents the UV cut-off for the effective theory. 

The approach of dimensional reduction is particularly valuable in disentangling contribu- 
tions from different degrees of freedom by accurate mixing analyses, as well as for treating larger 
temperatures T ^ which cannot be reached by 4d lattices. This can be used to inspect to 
what extent the plasma behaves perturbatively. fig. [6] then tells us that for any experimentally 
accessible temperature the largest correlation length of gauge-invariant operators belongs to 
the Aq gT degrees of freedom and not to the Ai g^T, in contrast to the naive parametric 
ordering eq. ()16|) . Only asymptotically is the perturbative ordering restored. Thus, on the soft 
scales gT and ~ g^T the plasma is strongly coupled. 



3.4 The free energy of static quarks 



The free energy of a static quark-antiquark pair is of interest for the physics of heavy quarkonia 
in the medium, and in particular to the question of J/T/i-suppression due to screening of the 
confining potential [47] . Potential models are then used to study bound state solutions, suitably 
generalized to finite temperature. More refined treatments try to establish a connection between 
the static potential and quarkonium spectral functions, to be discussed in the next section. Non- 
perturbative free energies serve as input for such applications, for recent work and references, 
see gH]. 

The qq free energy is defined [14] by the partition function of the thermal heat bath con- 
taining a static quark and antiquark source at separation r, 

( TrL(r) TtL\0) ) = cxp{-(F,, (|r|,T) - Fo(T) )/T}, (23) 

where Fq denotes the free energy of the heat bath. At zero temperature, i.e. Nt — > oo, Fqq 
reduces to the static quark potential. 

Let us begin our discussion in the pure gauge limit, where the expectation value of the 
Polyakov loop is a true order parameter. Numerical results below and above the pure gauge 
phase transition are given in [49] . When the temperature is switched on and increased towards 
Tc, a linearly rising free energy is maintained, implying infinite energy cost in separating the 
quark infinitely from the anti-quark, and thus confinement. However, the slope decreases with 
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Fig. 7. Left: Screening mass (quenched) for T > Tc from fits to eq. ([25}, d = 3/2 [49]. Right: Static 
quark free energy for Nf = 3 at temperatures 0.58 < T/Tc < 1.15 [53] ■ F{r) is normalized at (r = 1/T) 
to the T = Cornell potential, V{r)/^ = -a/r^/a + r^/a with a = 0.25 ± 0.05 (solid band). 



temperature, and the corresponding effective string tension is well fitted by the form 

a(0) 




(24) 



(predicted by an effective Nambu-Goto string (S^), with a = 1.21(5) and b = 0.990(5) ^49]. 
Above Tc an exponentially screened free energy is obtained, fitting an ansatz of the form 

In perturbation theory, the leading term originates from two-gluon exchange and predicts d — 2 
and exponential decay with twice the Debye mass [5T| . Lattice investigations [49l52j have shown 
that this simple behavior does not apply in the temperature and distance range explored. 
Rather, fits [49] to eq. ([25]) favoured d ~ 3/2 and found screening masses /i(T), shown in fig. [7] 
(left), to be compatible with the lowest color singlet 0^^ screening mass shown in figs. [5] [6] This 
is not surprising: the Polyakov loop is a gauge invariant operator, and since it is an exponential 
of gauge fields it couples to all J^'^ sectors. Consequently, its correlator decays with the lightest 
screening mass of the spectrum. 

When dynamical light quarks are present, the color charges of the heavy quarks are screened 
also below Tc and one observes [53l54l55j the expected string breaking, fig. [7] (right). The dis- 
tances where the free energies become flat in r range from 1.5 to 1 fm, decreasing with temper- 
ature, even at (bare) quark masses as large as m/T = 0.4. Note that the deviations from the 
zero temperature quenched potential set in already at distances of 0(0. 5fm) for temperatures 
> 0.75Tc. When the free energy is normalized to the short distance zero temperature potential, 
its large r asymptotic value rapidly decreases with T. 

The gauge invariant Polyakov loop averages over color singlet and octet configurations of 
the static source. Correspondingly, its correlator is often written as a superposition 14,,51..56j 

g-F,5(r,T)/T ^ 1 g-Fi(r,T)/T _,_ 8 ^_Fs(r,T)/T^ ^-26) 

with singlet and octet channels defined as 

e-^^(-^)/^ = i(Trit(x)L(y)), 

e-^^<^^'^y^ = i(Trit(x)TrL(y)) - l(TrLt(x)L(y)). (27) 
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While the singlet correlator is obviously gauge-dependent, the energy eigenvalues governing its 
decay in a spectral analysis are not ^57j and one might hope to gain useful insight into the 
colour dynamics to be used in potential models. However, this is dangerous. A careful spectral 
expansion making use of projection operators on the representations reveals that both, the 
singlet and the octet channel, decay as ^ ^„ |c„[C/]p exp{—En/T) [58], where the energies En 
represent the same static potential and its excitations which also contribute to the "average" 
free energy, while the matrix elements Cn[U] depend on the gauge or the operator used to 
describe the spatial string. Thus, the gauge-invariant spectral information is the same for all 
channels, while Fi, and their difference from Fqq are properties of the operators used, leaving 
their physical meaning in doubt. 



3.5 Dynamical quantities: quarkonium spectral functions and transport coefficients 

Quarkonium physics in the plasma can also be formulated more rigorously in a quantum field 
theoretical way. The physical excitations of the plasma with a given set of quantum numbers 
are encoded in the retarded Green's functions, or real time correlators, of operators carrying 
those quantum numbers [3] . In general the poles of such (momentum space) Green's functions 
will be complex, thus they are not eigenvalues of the Hamiltonian. If the imaginary part is 
small compared to the real part, one speaks of a quasi-particle excitation, with the real part 
interpreted as its mass and the imaginary part as its decay width, due to Landau damping in 
the plasma. 

Unfortunately, Monte Carlo methods only work for euclidean actions, and direct numer- 
ical simulations of real time correlators are impossible. However, one may try to statisti- 
cally exploit the information encoded in euclidean Green's function as follows. Let Gr{T) = 
X^x y t(^(''^>*)^V'(x7 i)V'(yi i + ''')rtp{y,t + t)) be a correlator of some meson operator in eu- 
clidean time. Then its Fourier transform is given by 

/•OO T 

Gr(r,p)=y^ ^Pr{cu,p)K{T,co), K{T,iu) = e^^nniuj) + e-^-[l + nsico)] , (28) 

where all the time dependence resides in the kernel K{t, uj), which contains only thermal distri- 
bution functions. All the dynamical, spectral information of the theory is contained in p(lu, p), 
which will not change under analytic continuation to Minkowski time. Still the difficulties are 
formidable. Computing Gr by lattice simulations, the left hand side is given only by a finite 
number of points and the inversion of the integral is an ill-defined problem. The way out is 
to "guess" with statistical methods the most likely p{uj,p) which will fit the formula. This is 
known as maximum entropy method (MEM), whose details are explained in [59]. Clearly, one 
needs as many points in the time direction as possible, hence anisotropic lattices with at ^ as 
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are useful. Another problem is that the MEM needs a prior, i.e. a model function has to be 
provided to constrain the fitted function. Usually one models spectral functions perturbatively 
at large lu. For a discussion of this and other systematic difficulties, see [60 . 

Early results based on this approach, obtained with quenched simulations, indicated that 
J /tp does not melt up to temperatures of 1.5 — 2Tc |61j . thus providing another reason to term 
the QGP 'strongly coupled'. There are now also dynamical simulations confirming this picture, 
as shown in fig. [8] (left) for Nf = 2 [62]. The temperatures parametrized by the shown NtS 
correspond to Tc — 2Tc. At low temperatures, the position of the first peak corresponds to 
the zero temperature mass of the J/tp. As temperature increases, the peak shrinks but keeps 
existing until about 1.7Tc, when it finally disappears. Hence the conclusion that "J/4' doesn't 
melt" until then (for another interpretation see 48J). No physical meaning can be attributed to 
the second peak, which is interpreted as a lattice-distorted part of the free particle continuum. 
This is illustrated in fig. [8] (right), where the spectral funtion for free particles in the scalar and 
vector channel are plotted in the continuum and on the lattice. The finite cut-off is faking a 
peak, which should disappear in the continuum limit. 

Other quantities of tremendous phenomenological interest are transport coefficients. In par- 
ticular, the shear viscosity is defined as the slope of another spectral function, 



Id 



0=0 



C(r) = / d^x{T:ki{T,^)7:ki{0,0)) = / ^ K{t,lo)p^^{t), (29) 



where nki is the traceless part of the spatial energy- momentum tensor. The numerical procedure 
to compute this on a lattice is the same as above. Here accurate information, and hence high 
resolution, on the low frequency tail is required. An early numerical attempt [64j as well as a 
recent more refined analysis [65' for the pure gauge theory give small viscosity to entropy ratios, 
77/s ~ 0.1 — 0.2, consistent with observations in the quark gluon plasma. However, it has been 
pointed out that extracting the low frequency part is firstly sensitive to modelling the spectral 
function [60j . and second intrinsically unstable, though this problem could be removed by a 
simple rescaling trick |66j . Building on this, ref. |65j suggested a new method to systematically 
improve on the spectral function by expanding it in terms of an orthogonal function system and 
fitting its coefficients. Hence, MEM calculations might soon offer some control over sytematics. 

Finally, a different approach to compute spectral functions is currently emerging, gener- 
alizing heavy quark effective field theory methods to real time correlation functions at finite 
temperature [67j . In the heavy quark limit a certain meson correlator reduces to a Wilson 
loop in Minkowski time. Thus, the static potential reappears in a real time framework, but 
it develops an imaginary part which is responsible for the damping and melting of the bound 
state [65. It has been shown in HTL perturbation theory that indeed a broadening spectral 
function is obtained in this approach [68j . The imaginary part of the potential has recently 
been confirmed non-perturbatively by classical lattice simulations |69j . It will be interesting to 
see if this approach can be generalized to the full quantum theory in the future, thus providing 
a bridge between field theory and potential models. 



4 The nature of the QCD phase transition 

In the previous sections we discussed the properties of QCD matter on either side of the quark 
hadron transition, but haven't yet addressed what the nature of that transition is. To obtain 
an unambiguous answer to this question is in fact a very difficult task. In statistical mechan- 
ics, phase transitions are defined as singularities, or non-analyticities, in the free energy as a 
function of its thermodynamic parameters. However, on finite volumes, free energies are always 
analytic functions and it can be proved that singularities only develop in the thermodynamic 
limit of infinitely many particles, or 1/ — > 00 [7D] . This is particularly obvious in the case of 
lattice QCD, whose partition function is a functional integral over a compact group with a 
bounded exponential as an integrand. It is thus a perfectly analytic function of T, /x, V for any 
finite V. Hence, a theoretical establishment of a true phase transition requires finite size scaling 
(FSS) studies on a series of increasing and sufficiently large volumes in order to extrapolate to 
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Fig. 9. Left: Conjectured phase diagram for Nf — 2 QCD with finite light quark massees. The physical 
case Nf = 2 + 1 is believed to qualitatively look the same, according to universality and continuity 
arguments as well as results from QCD-like models [2^. Right: Schematic phase transition behaviour of 
Nf =2 + 1 flavor QCD at fi = for different choices of quark masses {mu,d, nia), at fi = 0. 

the thermodynamic limit. Three different situations can emerge: a first order phase transition 
is characterized by coexistence of two phases, and hence a discontinuous jump of the order 
parameter (and other quantities) when going through such a transition. A second order transi- 
tion shows a continuous transition of the order parameter accompanied by a divergence of the 
correlation length and some other quantities, hke the heat capacity. Finally, a marked change 
in the physical properties of a system may also occur without any non-analyticity of the free 
energy, in which case it is called an analytic crossover. A familiar system featuring all these 
possibilities is water, with a weakening first oder liquid-gas phase transition terminating in a 
critical endpoint with Z(2) universality, as well as a triple point where the first order liquid-gas 
and solid-liquid transitions meet. Similar structures are also conjectured to be present in the 
QCD phase diagram [2|71| . fig. [9] (left). In this section we focus on the order of the transition 
for /i = as a function of Nf and quark masses, before tackling 7^ in sec. \5\ 

4.1 Universality and finite size scaling 

A most fascinating phenomenon in physics is the "universality" exhibited by physical systems 
near a critical point of second order phase transition. It is due to the divergence of the correlation 
length, which implies that the entire system acts as a coherent collective. Hence, microscopic 
physics becomes unimportant, the collective behaviour is determined by global symmetries. 
With the divergence of the correlation length the characteristic length scale disappears from 
the problem, and thermodynamic observables in the critical region obey scale invariant power 
laws. For example, at a second order ferromagnetic phase transition, the order parameter mag- 
netization vanishes as M ^ \T — T'cP, while the specific heat and correlation length diverge, 
C ~ \T — Tc\~°' and ^ ~ \T — Tc\~'^ , respectively. The critical exponents a, (3, v and similar ones 
for other quantities are determined by the effective global symmetries at the critical point, and 
thus are the same for all systems within a universality class. The latter are labeled by spin mod- 
els, since those can be readily solved numerically. The power law behaviour for thermodynamic 
functions follows from the scaling form of the singular part of the free energy. 



where L denotes the box size, and r and h are the relevant scaling fields. In the case of the Ising 
model these are the reduced temperature and magnetic field. Unfortunately, for most systems, 
and particularly QCD, it is not obvious which global symmetries the system will exhibit at a 
critical point, those are a dynamically determined subset of the total symmetry of the theory. 
Moreover, since there is no true order parameter in the case of finite quark masses, fields and 
parameters of QCD map into those of the effective model in a non-trivial way. 



ULy^T,Ly'^h), 



(30) 
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, 10. Finite volume scaling behaviour of specific heat and chiral susceptibility on Nt — 4. For 0(4) 
0(2) behaviour, the data should fall on a horizontal line [76] . 



4.2 The nature of the QCD phase transition for Nf = 2 + 1 at ^ = 0: qualitative picture 

Before discussing numerical results, let us briefly outline the qualitative picture for the order 
of the phase transition, fig. [9] (right). As mentioned in sec. 12.51 for Nf = 3 in the limits of zero 
and infinite quark masses (lower left and upper right corners) , order parameters corresponding 
to the breaking of a symmetry can be defined, and one finds numerically at small and large 
quark masses that a first-order transition takes place at a finite temperature Tc- On the other 
hand, one observes an analytic crossover at intermediate quark masses. Hence, each corner 
must be surrounded by a region of first-order transition, bounded by a second-order line. The 
situation is less clear for the chiral limit of the two flavour theory in the upper left corner. If 
the transition is second order, then chiral symmetry SUl{2) x SU{2)fi ~ 0(4) puts it in the 
universality class of 3d 0(4) spin models. In this case there must be a tricritical strange quark 
mass m*''*'^, where the second order chiral transition ends and the first order region begins. The 
exponents at such a tricritical point would correspond to 3d mean field [75. On the other hand, 
a first order scenario for the chiral limit of Nf = 2 is a logical possibility as well. 



4.3 The chiral transition for Nf = 2: is it 0(4) or first order? 

On the lattice, 0(4) will effectively look like 0(2) as long as there are discretization effects 
[72j . Wilson fermions appear to see 0(4) scaling ^73^, while staggered actions are inconsistent 
with both 0(4) and 0(2) [74]. (The staggered strong coupling limit, however, does display 
0(2) scaling [75]). An attempt to tackle this question by means of a finite size scaling analysis 
with unprecedented lattice sizes was made in [T^. The work simulates x 4 lattices with 
i = 16 — 32, using standard staggered fermions with several quark masses, the smallest being 
rn/T < 0.055. In a critical region, the specific heat or the chiral susceptibility scale as 

Cy - Co ~ L-^'^fc (tL^'\ amLy^y r = 1 - T/T^ 

X^L^'^f^irL^'^^amLy-). (31) 

Here the non-singular part of the specific heat Co has been subtracted. The authors of [7S^ thus 
fix Hh to its critical value, and then choose L and m for a series of simulations such as to keep 
{amLy^) constant, reducing the scaling problem on one variable only. The infinite volume limit 
in this procedure thus corresponds to the chiral limit and allows to check whether the data are 
consistent with the predicted scaling behaviour. 

Fig.flOlshows simulation results from [TSj. Scaling as in eq. (|3T|) would imply the data points 
to fall on a horizontal line, which is clearly not the case. After many variations on the analysis, 
the authors conclude that their data are closer to first order behaviour than to 0(4)/0(2). 
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Fig. 11. Left: The Binder cumulant as a function or quark mass for = 3, Nt —4 and different 
lattice sizes [82]. Right: The chiral critical line in the bare quark mass plane at ^ = 0, Nt = 4. Nf = 3 
is indicated by the solid line. Also shown are the physical point according to [85] , and an extrapolation 
using mean field exponents to the chiral limit, giving ml^^'^ ~ 2.8T. The point marked by the arrow 
has rriTr/mp = 0.148(2), compared to the physical value 0.18. From [82] , 



A different conclusion was reached in |77| , in which xQCD was investigated numerically. This 
is a theory modified by an irrelevant term (i.e. one going to zero in the continuum limit) such 
as to allow simulations in the chiral limit. The authors compare their data with those obtained 
from an 0(2) spin model on moderate to small volumes and find them to be compatible. 
This might point at finite volume effects of current Nf = 2 QCD simulations. The scaling 
region of the chiral limit might be excessively small, which would require extraordinarily large 
volumes in order to see the correct scaling. Another possibility for the failure of simulations 
to unambiguously pin down the order of this transition is of course the failure of the fermion 
discretizations to reproduce chiral symmetry. Future studies on either finer lattices or with 
chiral fermion actions should help to clarify this issue. 

Another possibility is to study the strength of the ?7a(1) anomaly discussed in sec. 12.51 
which also enters the effective sigma model [16] for the critical region, ft has recently been non- 
perturbatively demonstrated that a strong anomaly is required for the chiral phase transition 
to be second order [78] . 



4.4 Nf ~ 2 + 1: 3d Ising universality and the critical line m^(m„ d) 

Next we move our attention to the second order boundary line separating the first order region 
from the crossover in fig. [9| (right), starting with degenerate quark masses, Nf = 3. In this case 
on iV( = 4 lattices the critical quark mass marking the boundary between the crossover and 
the first order region is large enough for simulations to be carried out, and it was possible to 
determine it quite accurately |79|80|81(82] . together with the Z(2) universality class of the 3d 
Ising model to which it belongs [73]. The most practical observable to determine the latter is 
the Binder cumulant, constructed from moments of fluctuations of the order parameter, 

54(m,M)= //.^f^Xl ' 5Vi^ = ^V'-W), (32) 

evaluated on the phase boundary (3 = /3c(w, ^). In the infinite volume limit, B4 — > 1 or 3 for 
a first order transition or crossover, respectively, whereas it approaches a value characteristic 
of the universality class at a critical point. For 3d Ising universality one has B4 — 1.604 [83]. 
Hence is a non-analytic step function, which gets smoothed out to an analytic curve on 
finite volumes, with a slope increasing with volume to gradually approach the step function. 

A numerical example is shown in fig. [Tl] (left). The slope increases with lattice size, and the 
data fit the finite volume scaling formula predicted by universality, 



B4{m, L) =bo + hL^/^iam -~ am"). 



(33) 
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Fig. 12. Finite size scaling of the chiral susceptibility at the physical point for two lattice spacingsA^'t = 
4 (left) and Nt — 6 (right). The peaks saturate at a finite height, consistent with an analytic crossover. 

One observes that all curves intersect at the critical i?4-value, and moreover v is consistent 
with its 3d Ising value v = 0.63. Hence we can read off the critical quark mass in lattice units, 
am" = 0.0263(3) or rnVT^ = 1.052(1) [82]. 

While the nature of the transition and its universality class are determined by long range 
fluctuations and thus should be insensitive to the cut-off effects on a coarse lattice, the critical 
quark mass is a quantity requiring renormalization and is tremendously sensitive to such effects. 
Indeed, the critical quark mass and the corresponding pion mass have been computed on iVt = 4 
lattices with a standard staggered and a p4-improved staggered action. The critical pion mass 
in the improved case is almost a factor of four smaller than the non- improved result '79'34]. 
Hence, the location of the boundary line in the phase diagram is still far from the continuum 
result. 

Starting from the critical quark mass for Nf = 3, the boundary line has also been mapped 
for the non-degenerate Nf = 2 + 1 theory, again with A^t = 4 [52] ■ The question now arises 
which flavours to couple to the quark chemical potential. Within QCD alone there are no 
flavour changing interactions, while the initial state in heavy ion collisions has net u, d-quark 
numbers only. Thus it is sensible to assign a chemical potential to the two degenerate light 
quarks only. In nature, where the weak interactions are switched on as well, the situation is 
somewhat more complicated than this [84,. The methodology then is the same as in the three 
flavour case: fix a strange quark mass arus and scan the Binder cumulant in the light quark 
mass amu,d for the corresponding critical point. Repeating for several strange quark masses 
produces the critical line am'^{amu,d) shown in fig. [TT] (right). The results are in qualitative 
agreement with expectations from fig. [9] (right). They are also consistent with the possible 
existence of a tricritical point {rriu^d ~ 0,TOs = ml"'^). Using its known, mean field exponents, 
the data favour a heavy ml"'' ^ 2.8Tc. Again large corrections, presumably towards a smaller 
value, are expected in the continuum limit. Also, Tc is found to vary little along the critical 
line, in accordance with model calculations employing effective chiral lagrangians [86j . 

4.5 The nature of the transition at the physical point 

A more immediate issue is whether the QCD physical point indeed lies on the crossover side 
of the critical line as expected. For that purpose spectrum calculations at T ~ O have been 
performed ^82] at the parameters indicated by the arrows in fig. [11] (right) . They show that TOs 
at the upper arrow is approximately tuned to its physical value (^^ ~ ^^|phys=0.65), while 
the pion is lighter than in physical QCD (^ = 0.148(2) < ^Iphys = O.lsj.^This confirms that 
the physical point lies to the right of the critical line, i.e. in the crossover region. 

An important question is whether this finding, calculated on coarse lattices Nt — 4, is 
stable under continuum extrapolations. This has been recently addressed and affirmed in [87j . 
using a Symanzik improved gauge and a stout-link improved staggered fermion action for four 
different lattice spacings, Nt — 4,6,8,10. The simulations were carried out along the lines of 
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constant physics, i.e. the bare parameters were tuned such that the hadron spectrum in physical 
units is constant for different lattice spacings. The result of a finite size scaling analysis of the 
susceptibility of the chiral condensate is shown in fig. [121 The peak height appears to clearly 
saturate at some finite value, and the transition thus corresponds to a crossover. The finite 
peak height was calculated for all four lattice spacings and extrapolates to a finite value also 
in the continuum limit. Hence, one concludes that the transition between hadronic and quark 
gluon plasma physics at zero density indeed is an analytic crossover for physical QCD. In view 
of the concerns about the validity of the staggered fermion actions |8j, it would of course be 
desirable to reprodce this result with alternative fermion discretizations. 



5 Lattice QCD at finite density 

Since 2001, significant progress was made towards simulating QCD with realistic parameter 
values at small baryon densities. This is achieved by a number of methods that circumvent 
the sign problem, rather than solving it: i) Multi-parameter reweighting, ii) Taylor expansion 
around ^ — and iii) simulations at imaginary chemical potential, either followed by analytic 
continuation or Fourier transformed to the canonical ensemble. It is important to realize that 
all of these approaches introduce some degree of approximation. However, their systematic 
errors are rather different, thus allowing for powerful crosschecks. All methods are found to be 
reliable as long as fJ,/T < 1, or /zs < 550 MeV, which includes the region of interest for heavy 
ion collisions. Reviews specialized on this subject can be found in |88I89| . 



5.1 Two parameter reweighting 

Significant progress enabling finite density simulations was made a few years ago, by a gen- 
eralization of the Glasgow method [90] to reweighting in two parameters [91]. The partition 
function is rewritten identically as 



where the ensemble average is now generated at /i = and a lattice gauge coupling f3o, while 
a reweighting factor takes us to the values of interest. The original Glasgow method |90] 
reweighted in /j, only and was suffering from the overlap problem: while the reweighting formula 
is exact, its Monte Carlo evaluation is not. The integral gets approximated by a finite number 
of the most dominant configurations, which are different for the reweighted and the original 
ensemble, and this difference grows with fi. When searching for a phase transition at some 
fj,, one-parameter reweighting uses a non-critical ensemble at /i = 0, thus missing important 
dynamics. By contrast, two-parameter reweighting uses an ensemble generated at the pseudo- 
critical coupling Pci^J■ = 0), which is then reweighted along the pseudo-critical line of the phase 
change. Thus one is working with an ensemble that probes both phases, improving the overlap 
with the physical ensemble. 

A difficulty in this approach is that the determinant needs to be evaluated exactly, which is 
costly. Moreover, because of the sign problem the reweighting factor is exponentially suppressed 
with volume and chemical potential, thus limiting the applicability to presently moderate values 
of those parameters. Moreover, since the statistical fluctuations are those of the simulated 
ensemble instead of the physical one, all reweighted measurements are correlated and it is 
difficult to obtain reliable error estimates. For a proposed procedure see [92]. A further problem 
arises with staggered quarks, where the root of the determinant has to be taken. For fJ, ^ this 
may enhance cut-off effects to 0{a) instead of O(a^) for /x = [93] . 





(34) 
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5.2 Finite density by Taylor expansion 

Another way to gain information about non-zero /i is to compute the coefficients of a Tay- 
lor series expansion of observables in powers of ^/T . Early attempts have looked at suscep- 
tibilities and the response of screening masses to chemical potential |94|95|96|97j . More re- 
cently it has also been used to gain information on the phase transition and its nature itself 
[98199110011011102] , This idea exploits the fact that on finite volumes the partition function 
Z{m > 0,^,T) is an analytic function of the parameters of the theory. For small enough fj,/T 
one may then hope to get away with only a few terms, whose coefficients are calculated at 
fi = Q. Moreover, because of CP symmetry of the QCD action the partition function is even 
in /i, Z(^) = Z(—^), such that physical observables have series expansions in (/i/T)^. Thus, 
e.g. the pressure density can be expressed as an even power series, 

p{T,^,) = -^=(pj\ogZ{T,^,)^ |_ = ge2„(T)(^)'". (35) 

^ ^ n=0 

The coefficients are equivalent to generalized quark number susceptibilities at /Lt = and measur- 
able with standard simulation techniques. Since all the /i-dependence of the partition function 
is in the fermion determinant, its derivatives need to be computed, 

tr (^Af-i^A/-i^ , etc., (36) 

and iterations for higher orders. These expressions become increasingly complex and methods 
to automatize their generation have been devised I102j . Note that one now is dealing with traces 
of composite local operators, which greatly facilitates their numerical evaluation by statistical 
estimators compared to a computation of the full determinant as required for reweighting. 

For high enough temperatures T>Tc, the scale of the finite temperature problem is set by 
the Matsubara mode ~ ttT, and one would expect coefficients of order one for an expansion 
in the 'natural' parameter ^/{ttT) [81\. This is indeed borne out by numerical simulations and 
explains the good convergence properties observed for /i/T< 1. 



aindetM _ (m-^^^'^\ 



5.3 QCD at imaginary ^ 

The hermiticity relation eq. tells us that the QCD fermion determinant with imaginary 
fi — i^i is real positive, hence it can be simulated just as for = 0. It is then natural to 
ask whether such simulations can be exploited to learn something about physics at real /z. 
The strategy to get back to real /i is to fit the full Monte Carlo results to a truncated Taylor 
series in fii/T. In case of apparent convergence it is then easy to analytically continue the 
power series to real fj,. This idea was first used for observables like the chiral condensate and 
screening masses in the deconfined phase |103I96| . It was then shown to be applicable to the 
phase transition itself JL04J i which has recently been exploited in a growing number of works 
j81l82ll05ll06ll()7ll08ll()9] . 

For complex /i = Hr + ifJ^i, the QCD partition function eq. ([T]) is periodic in the imagi- 
nary direction, with period 2tt/Nc for Nc colours [llOj . Hence, in addition to being even in fx, 
the partition function has the additional symmetry Z(fir/T, fXi/T) = Z{fj,r/T^ fii/T + 27r/3). 
Because of the anti-periodic boundary conditions on fermions, this symmetry implies that a 
shift in fii by certain critical values is equivalent to a transformation by the Z{3) centre of 
the gauge group. Thus, there are Z{3) transitions between neighbouring centre sectors for all 
{pLi/T)c = ^ (n -I- i) , n = 0, ±1, ±2, .... It has been numerically verified that these transitions 
are first order for high temperatures and a smooth crossover for low temperatures |104|105| . 
This limits the radius of convergence for analytic continuation, which is given by the first such 
transition, or ji/T = 7r/3. Hence the approach is limited to |/Lt/T| < 1. 
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Within this circle, the strategy then is to measure expectation values of observables at 
imaginary fi and fit them by a Taylor series, 



Working at imaginary /i has a couple of technical advantages. It is computationally simple and 
much cheaper than reweighting or computing coefficients of the Taylor expansion. Moreover, 
both parameters /?, /i are varied and thus one obtains information from statistically indepen- 
dent ensembles. It also offers some control on the systematics by allowing a judgement on the 
convergence of the fits. Furthermore, it is a good testing ground for effective QCD models: 
analytic results can always be continued to imaginary /i and be compared with the numerics 
there, as demonstrated for several examples in [106j . 

A different approach making use of imaginary chemical potential is to employ the canonical 
ensemble at fixed quark number, which is related to the grand canonical ensemble via the 
integral transform 



One may then compute the grand canonical partition function at imaginary /i and perform the 
Fourier transform numerically ^lllj . This will only work for moderate Q and small volumes, 
as the sign problem is reintroduced by the oscillations of the exponential. For small lattices 
6^ X 4, 4^, first interesting results on the phase diagram have been obtained, both for staggered 
[112| and Wilson fermions |113| . However, this approach has no overlap problem, and one might 
hope to push to larger chemical potentials once computational resources are available. 



5.4 Plasma properties at finite density 

Having developed computational tools for finite density, one can repeat the studies discussed in 
the previous sections and see how finite baryon densities affect the screening masses [34195196197] , 
the equation of state |100|101fTT4j or the static potential 115j. In all those cases the influence 
of the chemical potential is found to be rather weak, and for lack of space we will not further 
discuss those calculations here. We likewise have to pass over the work done on certain modifi- 
cations of QCD, for which the sign problem is manageable or absent, such as QCD in the static 
limit, two-colour QCD or QCD at finite isospin. Such studies give important qualitative insights 
and some are covered in previous reviews |116|88j . Instead we concentrate here on calculations 
of the QCD phase diagram at finite density, where the order of the phase transition is expected 
to change as /i is increased, fig. [9] (left). 



5.5 The critical temperature at finite density 

As in the case of zero density, let us first discuss the phase boundary, i.e. the (pseudo-) critical 
temperature Tc{fJ,), before dealing with the order of the phase transition. The critical line 
has been calculated for a variety of flavours and quark masses using different methods. For 
a quantitative comparison one needs data at one fixed parameter set and also eliminate the 
uncertainties of setting the scale. Such a comparison is shown for the critical coupling in fig. 1131 
(left), for Nf = 4 staggered quarks with the same action and quark mass m/T sa 0.2. (For 
that quark mass the transition is first order along the entire curve) . One observes quantitative 
agreement up to /i/T « 1.3, after which the different results start to scatter. This vindicates our 
earlier statement that all methods appear to be reliable for n/T < 1. Note that strong coupling 
results at /? = predict a^dP = 0) = 0.35 — 0.43 [117], requiring the line to bend down rapidly, 
and thus favour the data from the canonical approach. 
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Fig. 13. Left: Comparison of different methods to compute the critical couplings [112j . Right: The 
phase diagram for physical quark masses as predicted by the two parameter reweighting method [85] . 
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Table 1. Coefficient t2 in the Taylor expansion of the transition line, eq. (|39p . All results have been 
obtained with Nt = 4. 



The case of physical quark masses, after conversion to continuum units, is shown in fig. 1131 
(right) [85] . One observes that the critical temperature is decreasing only very slowly with /i. 
This is consistent with a description by a Taylor expansion in powers of (^/ttT)^ with coefficients 
of order one, 

m-^-<^i^,-'nM^r+o((JLY] . (39) 



Tc{0) ' ^'\ttTJ \\ttTJ 

The leading coefficients for various cases have been collected from the literature in [55] and 
are reproduced in table [1] The curvature gets stronger with increasing Nf, which is consistent 
with the observation at zero density that T^, is lowered by light flavours, cf. sec. 13.11 Subleading 
coefficients are also beginning to emerge but at present not statistically significant yet. It has 
been noted that the convergence should speed up when constructing Pade approximants for 
the series, which also tends to increase the curvature towards larger /i [1061118] . Finally one 
should note that the continuum conversions relying on the two-loop beta function are certainly 
not reliable for these coarse lattices, while fits to non-perturbative beta functions tend to also 
increase the curvature. 



5.6 The QCD phase diagram for /i 7^ and the critical point 

As in the case of /i = 0, a determination of the order of the transition, an hence the search 
for the critical endpoint, is much harder, and we begin by discussing the qualitative picture. 
If a chemical potential is switched on for the light quarks, there is an additional parameter 
requiring an additional axis for our phase diagram characterizing the order of the transition, 
fig. [9] (right). This is shown in fig. [M] where the horizontal plane is spanned by the /i = 
phase diagram in ms,mu^d and the vertical axis represents /z. The critical line separating the 
first order section from the crossover will now extend to finite fi and span a surface. A priori it 
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Fig. 14. Upper panel: The chiral critical surface in the case of positive (left) and negative (right) 
curvature. If the physical point is in the crossover region for ^ = 0, a finite phase transition will only 
arise in the scenario (left) with positive curvature, where the first-order region expands with [s,. Note 
that for heavy quarks, the first-order region shrinks with [i, cf. sec. 15.91 Lower panel: phase diagrams 
for fixed quark mass (here Nj — 3) corresponding to the two scenarios depicted above. 



is, of course, not known whether this surface might be leaning to the right or the left, or even 
have a more complicated behaviour as a function of the quark masses. However, the expected 
QCD phase diagram is only obtained if the left situation of fig. [T3]is realized (unless there are 
additional critical surfaces yet unknown). The first order region expands as fi is turned on, 
so that the physical point, initially in the crossover region, eventually belongs to the critical 
surface. At that chemical potential fiE^ the transition is second order: that is the QCD critical 
point. Increasing fi further makes the transition first order. A completely different scenario 
arises if instead the first-order region shrinks as fi is turned on. In that case the physical point 
remains in the crossover region for any /z, fig. [14] (right). 

There are then different strategies to determine the QCD phase diagram. One can fix a 
particular set of quark masses and for that theory switch on and increase the chemical potential 
to see whether a critical surface is crossed or not. Such calculations are covered in sec. 15.71 
Alternatively, sec. 15.81 discusses how to start from the known critical line at fi = and study 
its evolution with a finite fi. That will give information on the whole phase diagram in fig. I14[ 
including, eventually, physical QCD. 



5.7 Critical point for fixed masses from reweighting and Taylor expansion 

Reweighting methods produced the first finite density phase diagram from the lattice, ob- 
tained for light quarks corresponding to m.^ ^ 300 MeV jll9j . A later simulation at physical 
quark masses puts the critical point at fig ~ 360 MeV [85], fig. [15] (left). In this work x 4 
lattices with L = 6—12 were used, working with the standard staggered fermion action. 
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Fig. 15. Left: Lee- Yang zeroes in the complex /3-plane for SU(3) pure gauge theory |120| . Right: 
Imaginary part of the Lee- Yang zero closest to the real axis as a function of chemical potential [85) . 



Quark masses were tuned to m„^d/7c ~ 0.037, rris/Tc ~ 1, corresponding to the mass ratios 
rriTr/mp « 0.19, mT^/mK ~ 0.27, which are close to their physical values. A Lee- Yang zero anal- 
ysis |70| was employed in order to find the change from crossover behaviour at /i = to a first 
order transition for fi > fj,E- This is shown in fig. 1151 For a crossover the partition function has 
zeroes only off the real axis, whereas for a phase transition the zero moves to the real axis when 
extrapolated to infinite volume. For a critical discussion of the use of Lee- Yang zeros in combi- 
nation with reweighting, see [120j . Recently, reweighting has been combined with the density of 
states method [121| . in order to extend the applicability of reweighting to larger values of fi/T. 
First interesting results have been obtained for Nj — 4, indicating a new high density phase 
and a possible triple point, where the high density transition line meets the deconfinement line 
[121] . Unfortunately, the method is computationally very expensive and so far limited to coarse 
and small lattices, so that it is difhcult to unambiguously establish those findings at present. 

In principle the determination of a critical point is also possible via the Taylor expansion. 
In this case true phase transitions will be signalled by an emerging non-analyticity, or a finite 
radius of convergence for the pressure series about /j, = as the volume is increased. The 
radius of convergence of a power series gives the distance between the expansion point and 
the nearest singularity, and may be extracted from the high order behaviour of the series. A 
possible definition is 

l/2n 

p = lim pn with pn = 



C2n 



(40) 



General theorems ensure that if the limit exists and asymptotically all coefficients of the series 
are positive, then there is a singularity on the real axis. More details as well as previous 
applications to strong coupling expansions in various spin models can be found in [122j . In 
the series for the pressure such a singularity would correspond to the critical point in the 
(/i, r)-plane. 

A critical endpoint for the Nf = 2 theory, based on this approach, was reported in |102j . 
The authors perfomed simulations on x 4 lattices with L — % ~ 24, using the standard 
staggered action. The quark mass was fixed in physical units to m/Tc = 0.1. The aim of the 
simulations was to bracket the critical point by computing the Taylor coefficients of the quark 
number susceptibility up to sixth order (i.e. 8th order for the pressure) for various temperatures 
in the range T /Tc = 0.75 — 2.15, and extrapolate to finite p. This was done for different lattice 
volumes in order to get an estimate of finite voulme effects. 

The results for the convergence radius eq. (|40p are shown in fig. [16] (left). A rather strong 
volume dependence is apparent. While for the smaller 8^ lattice the estimator p„ does not 
seem to converge to a finite radius of convergence, the results on the larger 24^ lattice are 
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Fig. 16. Left: Estimators of the radius of convergence p„, eq. (|40p . at T/Tc = 0.95 on A'^t = 4 lattices. 
Circles represent L = 8, squares L = 24 [102] . Right: Quark number susceptibility computed through 
0(/i'') (dashed lines) and 0(/i'') (solid lines) [TOT] . 



consistent with settling at a limiting value. Taking the large volume result at face value and 
extrapolating to all orders the estimate for the location of the critical point is /i^/Tg = 1.1 ±0.2 
at Te/Tc{h = 0) = 0.95 \W2\ . 

Another investigation of the two-flavour theory, also using the Taylor expansion of the 
pressure, is reported in [1001101 ]. This group works with a 16"^ x 4 lattice with p4- improved 
staggered fermions and a Symanzik-improved Wilson action, the quark mass is set to m/T w 0.4. 
The calculation to order /Lt* was performed in [lOOj while results on /i^ are presented in [lOlj . The 
last work also contains detailed discussions of analytic calculations to compare with, namely 
the pressure in high temperature perturbation theory [123 , which is going to hold at very high 
temperatures, as well as the hadron resonance gas model, which gives a rather good description 
of the pressure in the confined phase [124j . 

In agreement with [102j and qualitative expectations, their detailed results for the coeffi- 
cients in the pressure series satisfy cg <C C4 <C C2 for T > Tc, i.e. one would have coefficients 
of order one for an expansion in (fi/irT). An impression of the convergence of the series can 
be obtained by looking at the quark number susceptibility calculated to consecutive orders, as 
shown in fig. [16] (right). For T< 1.2Tc, the series seems to converge rapidly and the /i^-result 
is compatible with the one through order /z^. Around the transition temperature Tc, the /i"*- 
results show a peak emerging with growing /i/Tc, which in [100] was interpreted as evidence 
for a critical point. However, the 11^ contribution suggests that in this region results do not yet 
converge, and the structure is hence not a significant feature of the full pressure. Furthermore, 
estimates for the radius of convergence through that order agreed with predictions from the 
hadron gas model, which however has infinite radius of convergence. Hence the conclusion in 
[101| that there is no signal for a critical point at that quark mass. 



5.8 The change of the critical line with /i 

Rather than fixing one set of masses and considering the effects of ^, one may map out the 
critical surface in fig. [T3]by measuring how the fj, = critical boundary line changes under the 
influence of fi. This question was adressed using an imaginary chemical potential in [81|82j . The 
methodology employed is the same as in sec. 14.41 i.e. measurement of the Binder cumulant of 
the chiral condensate. The results for Nf = 3 are summarized in fig. [T7] (left). The chemical 
potential is found to have almost no influence on B4 as a function of quark mass. A lowest-order 
fit, linear in am and (a/i)^, gives the error band fig. 1171 (right), corresponding to 

am''{a^l) = 0.0270(5) - 0.0024(160) (a^)^ . (41) 

Care must be taken for the conversion to physical units. The crucial point is that the gauge 
coupling (3 is tuned with changing to stay on the critical line, so that a(/3) decreases. Ex- 
pressing the change of the critical quark mass with chemical potential in lattice and continuum 
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Fig. 17. Left: B4{am,a^i) for Nf — 3, Nt = i and imaginary chemical potentials. Right: One-sigma 
error band for the critical mass am'^^ayLi) resulting from a linear fit to the data on the left. From [82j . 



units as 



then ci and c'j^ are related by 

Since c'j^ is observed to be nearly zero, it is the second term that dominates, leading to an overall 
negative coefficient ci = —0.7(4) [82] , 

This is evidence that, in the Nf — 3 theory on an A^t = 4 lattice, the region of first- 
order transitions shrinks as a baryon chemical potential is turned on, and the "exotic scenario" 
of fig. [14] (right) is realized. Interestingly, similar qualitative conclusions are obtained from 
simulations of the three flavour theory with an isospin chemical potential [125] . as well as 
simulations at imaginary fi employing Wilson fermions |109j . 

This investigation has also been extended to the case of non-degenerate quarks [S^. Fig.fTSl 
(left) shows a comparison of the critical line for /i = with some critical masses calculated for 
a sizeable baryon chemical potential. The data show the same trend as for Nf = 3: the critical 
mass is nearly unchanged or slightly increasing with /i^, in lattice units. The conversion to 
physical units proceeds as in the three flavour case and gives a dominant negative contribution 
to ci . Together with a very small value for c'l , it implies again that the first-order region shrinks 
as the baryon chemical potential is turned on, and the "exotic scenario" of fig. [U] (right) is the 
correct one. 



5.9 The heavy quark limit: Potts model 

In light of these surprising results, it is also interesting to consider the heavy quark corner of 
the schematic phase diagram of fig. 1141 Simulations of quenched QCD at finite baryon number 
have been done in 126J. As the quark mass goes to infinity, quarks can be integrated out and 
QCD reduces to a gauge theory of Polyakov lines. First simulations of this theory with Wilson 
valence quarks can be found in |127j . At a second order phase transition, universality allows 
us to neglect the details of gauge degrees of freedom, so the theory reduces to the 3d three- 
state Potts model, which is in the 3d Ising universality class. Hence, studying the three-state 
Potts model should teach us about the behaviour of QCD in the neighbourhood of the critical 
line separating the quenched first order region from the crossover region. For large fx the sign 
problem in this theory has actually been solved by means of cluster algorithms [128j . 
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Fig. 18. Left: Comparison of the critical line at /i = and a/xj = 0.2. Right: The critical heavy quark 
mass separating first order from crossover as a function of fj,^ from the Potts model [129] . 



Here we are interested in simulations at small fi/T |129| . In this case, the sign problem is 
mild enough for brute force simulations at real /i to be feasible. In [129j . the change of the 
critical heavy quark mass is determined as a function of real as well as imaginary /x, as shown 
in fig. [TH] (right). Note that M'^{fj,) rises with real chemical potential, i.e. the first order region 
in fig. [14] shrinks as finite baryon density is switched on. This system is thus an example of the 
non-standard scenario discussed in the previous section. The calculation also gives some insight 
in the problem of analytic continuation: while fig. [15] clearly endorses the approach in principle, 
an 0(/i^)-fit was required to reproduce the data on both sides of /i^ — 0. Similar accuracy is 
much more difficult to achieve around the chiral line with present resources. 



5.10 Discussion of critical end point results 



The results about the critical surface from sec. 15.81 appear to be in qualitative contradiction 
with those of [85], |102j . which both conclude for the existence of a critical point {ij,e,Te) at 
small chemical potential ^e/Te 1- However, in considering the reasons for such disagreement, 
one can see that the different data sets are actually not inconsistent with each other, and the 
differing pictures can be explained by standard systematic effects. 

In '102] the critical point was inferred from an estimate of the radius of convergence of the 
Taylor expansion of the free energy. Regardless of the systematics when only 4 Taylor coefficients 
are available, the estimate is made for Nf — 2. The (/i, T) phase diagram of this theory might 
well be qualitatively different from that of Nf = 2 + 1 QCD, as illustrated in fig. [T9| (right). 
Its critical endpoint point, obtained by intersecting a critical surface when going up vertically 
from the Nf = 2 quark mass values, is clearly a long way from the critical endpoint of physical 
QCD, and nothing follows quantitatively from the value of one for the other. 

In |5S] the double reweighting approach was followed. By construction, this reweighting is 
performed at a quark mass fixed in lattice units: amu,d = ^-^^^ ~ const. Since the critical 
temperature Tc decreases with fi, so does the quark mass. This decrease of the quark mass 
pushes the transition towards first order, which might be the reason why a critical point is 
found at small fi. This effect is illustrated in the sketch fig. [T9| (left), where the bent trajectory 
intersects the critical surface, while the vertical line of constant physics does not. Put another 
way, [85] measures the analogue of c'^ instead of ci in eq. ([42]) . From fig. [T5l (right) we see that 
the distance of the physical theory from criticality stays initially constant, consistent with a 
coefficient c'l « as that in [82], sec. 15.81 Taking the variation a(//) into account could then 
make a dominant contribution, which might possibly change the results qualitatively. 

Conversely, in ^82j only the leading order coefBcient of to'^(/u) has been determined from 
imaginary ^. It cannot yet be excluded that there are cancelling terms of higher order, alter- 
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Fig. 19. Left: Effect of keeping the quark mass fixed in lattice units while changing fi, as in [85) . Right: 
Location of critical points for the Nf =2 + 1 and the Nf = 2 theory, as considered in [102'. Knowing 
the latter does not predict the former. 



nating in sign. Such contributions would no longer cancel after continuation to real /i, leading 
to a different picture. Konwledge of the next term in the series will clarify this. 

On the other hand, a robust finding is the high quark mass sensitivity of the critical point: 
irrespective of the sign, if ci ^ 0{1) in eq. ([42]), m'^{p) is a slowly varying function of just as 
the pressure, screening masses or T^- Hence /Z£;(to) is rapidly varying. A change of quark masses 
by a few percent will then imply a change of he by 0(100%). One should also remember that 
most investigations so far have used unimproved staggered quarks on coarse Nt — A lattices 
only. This might be worrisome given the exceedingly light quarks involved, cf. sec. 12.21 and 
[8193] . Finally, a more complicated picture with additional critical surfaces in fig. [T3] is yet 
another possibility. In light of these circumstances even the qualitative features of the QCD 
phase diagram cannot be regarded as settled until they are probed with better accuracy on 
finer lattices. 



6 Conclusions 

We have summarized our current understanding of finite temperature and density QCD from 
numerical studies on the lattice. While many results in the pure gauge theory are available in 
the continuum limit, simulations with dynamical fermions still suffer from systematic errors. 
These are mainly due to the finite lattice spacing, and quark masses which don't meet their 
physical values yet. For the critical temperature and the equation of state these problems are 
now being tackled and the first continuum extrapolations become available. 

Existing results provide us with a detailed picture of how equilibrium plasma properties 
change through the deconfinement transition up to a few Tc. Combinations of perturbative 
calculations and numerical simulations have produced insight into the regime of very high 
temperatures as well as the dynamics and mixing of hard and soft momentum modes. Altogether 
this provides a quantitative understanding of the relevant length scales in the plasma, as well 
as tests of the applicability of thermal perturbation theory. 

The naive picture of the deconfined phase as a weakly interacting parton gas is not sup- 
ported. For temperatures relevant to heavy ion collisions, the plasma displays strong residual 
interactions through soft gluonic modes, which cannot be treated perturbatively, and which in- 
fiuence different quantities in different ways. This gives a consistent explanation to the various 
observed features: the equation of state, susceptibilities and fermionic correlators are dominated 
by hard modes and not far from their ideal gas values, but the corrections are significant. Glu- 
onic correlators, on the other hand, are dominated by soft modes and entirely off their leading 
perturbative predictions. An ideal gas is established only at asymptotically high temperatures. 

Significant progress was made over the last five years regarding dynamical quantities as well 
as finite density simulations. In both fields methods have been developed, that are currently 
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being scrutinized for their reliability. Calculations of spectral functions provide a picture of 
remnant binding forces in the plasma phase as well as first semi-quantitative results for the 
transport coefficients. Calculations at finite density are now possible for /x/T<l, with good 
agreement between all methods whenever equal parameter sets arc compared. However, the way 
to a quantitative understanding of the QCD phase diagram is still long, mainly due to the high 
quark mass and cut-off sensitivity of the critical endpoint. Current developments, also regarding 
computing speed-up for different fermion discretizations and improvement programmes, give 
reason to believe that these questions can be significantly improved upon in the near future. 
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